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Spontaneous spin vortex formation in the magnetic phase transition of a trapped spin-1 Bose- 
Einstein condensate is investigated using mean-field theory. In a harmonic trapping potential, an 
inhomogeneous atomic density leads to spatial variations of the critical point, magnetization time 
scale, and spin correlation length. The Kibble-Zurek phenomena are shown to emerge even in such 
inhomogeneous systems, when the quench of the quadratic Zeeman energy is fast enough. For slow 
quench, the magnetized region gradually expands from the center of the trap pushing out spin 
vortices, which hinders the Kibble-Zurek mechanism from occurring. A harmonic trap with a plug 
potential is also taken into account. 



I. INTRODUCTION 

Symmetry breaking phase transitions are considered 
to play crucial roles in the early universe. As the hot 
universe cooled down, the phase transitions broke the 
symmetries of the vacuum fields. Since causally discon- 
nected regions acquire independent values of the order 
parameter in the course of the phase transition, topo- 
logical defects can be left behind 1 , such as monopoles, 
strings, and domain walls. It was proposed that this cos- 
mological scenario of topological defect formation can be 
tested by the normal fluid-supernuid phase transition of 
liquid helium 2 -. Such a mechanism of topological defect 
formation is called Kibble-Zurek (KZ) mechanism, which 
has been studied in a wide variety of systems 3 - - — . 

Bose-Einstein condensates (BECs) of atomic gases are 
highly controllable quantum systems and suitable for 
studying the KZ mechanism in a controlled manner. The 
BEC transition breaks the U(l) symmetry for a single- 
component system, and quantized vortices can be formed 
by the KZ mechanism. This has been demonstrated in 
the experiments reported in Refs j 12 i 13 . Spinor BECs 
(BECs of atoms with spin degrees of freedom) have a 
rich variety of magnetic phases with different symmetry 
groups, and thus have various kinds of topological de- 
fect o 14 i 15 . In the experiments reported in Refs j 16 i 17 , the 
transition from the polar state to the ferromagnetic state 
in a spin-1 87 Rb was observed, which was controlled by 
an external magnetic field. Formation of spin vortices by 
the KZ mechanism in this magnetic transition has been 
investigated in Refs^~— . It is predicted that the KZ 
mechanism can also be tested by the Mott transition of 
cold atoms in an optical lattice^ 2 -, soliton formation in the 
BEC transition in a one-dimensional gas 23 , a miscible- 
immiscible transition in a binary BEC 24 , and a magnetic 
transition in an antiferromagnetic spinor BEC 25 . 

In Ref. 21 , we studied the KZ mechanism in the mag- 
netic transition of a spin-1 87 Rb BEC and numerically 
demonstrated the KZ scaling properties. However, the 
numerical simulations in Ref *2I were restricted to the sys- 
tems with uniform atomic density. In the present paper, 
we perform numerical simulations of the magnetization 



dynamics of a spin-1 BEC confined in a harmonic trap- 
ping potential to show that the KZ mechanism can be 
observed in realistic experiments. We will show that the 
inhomogeneity of the trapped system has two effects on 
the KZ properties. The first one is caused by the spatial 
dependence of the spin correlation length. The number 
of spin vortices created by the KZ mechanism depends 
on the spin correlation length, and therefore, depends on 
the position. The second one originates from the compe- 
tition between two velocities. Since the density is high 
around the center of the atomic cloud, the magnetiza- 
tion starts from the center and the magnetized region ex- 
pands outward. If this expansion velocity is slower than 
the velocity of the spin wave, the magnetized region can 
be causally connected with the region that is going to 
magnetize, and the KZ mechanism breaks down. A plug 
potential applied to the center of the trap is shown to 
resolve this problem. 

This paper is organized as follows. Section [TT] formu- 
lates the problem and provides mean-field and Bogoli- 
ubov analyses. Sections IIII Al and IIII Bl show the numer- 
ical results for sudden quench and gradual quench of the 
magnetic field, respectively. Section IIII CI examines the 
case of a harmonic potential with a plug potential. Sec- 
tion |IV] concludes this paper. 



II. MEAN-FIELD ANALYSIS OF SPIN 
CORRELATIONS 

We consider bosonic atoms with mass M and hyperfine 
spin F = 1 confined in an external potential Vt ra p(r). 
The magnetic field B is applied in the z direction, and the 
linear and quadratic Zeeman effects change the energies 
of spin sublevels m = ±1 by 



P = TQfUbB, 
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(1) 



respectively, where gp is the hyperfine g factor, /iB is the 
Bohr magneton, and E^f is the hyperfine splitting energy. 
For 87 Rb atoms, g F = 1/2 for F = 1 and E M /h ~ 6.8 
GHz. The interaction between atoms is characterized by 
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spin-independent and spin-dependent interaction coeffi- 
cients given by 



c 



~M 3 : 
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respectively, where as is the s-wave scattering length for 
two colliding atoms with total spin S. We use the values 
of a = 101.8a B and a 2 = lOOAa^ for F = 1 87 Rb 
atoms, where <2b is the Bohr radius. 

We employ the mean- field theory at zero temperature. 
The state of the system is described by the macroscopic 
wave functions ^> m (r,t). The mean- field energy is given 

by 



trap + mp + m ^ m 
(3) 



where 



p(r,t) = l^p + l^oP + IV'-il 

^X*",*) = 51 i'mfmm'lpm', 
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with / = (/ x , / y , / z ) being the spin-1 matrices. The dy- 
namics is given by the Gross-Pitaevskii (GP) equation, 



ih 



dt 
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where the right-hand side indicates functional deriva- 
tive. In the rotating frame of the spin space (^±i 
e^ ipt / n il)±i), the linear Zeeman terms in Eq. ([6]) can be 
eliminated, and we neglect them in the following calcu- 
lations. 

To study the behaviors of the system analytically, we 
consider a uniform system with density p = uq. When 
ci < and q > 0, which is the case of spin-1 87 Rb, the 
ground state of Eq. (|3]) satisfying F z = is given by 
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for q > q c and 
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for q < q Cl where 



2|ci|n , 



(8) 
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and a and /? are arbitrary phases. The states in Eqs. ([7]) 
and ((SJ) are called the polar state and broken axisym- 
metry state 27 , respectively. The transverse magnetiza- 
tion of the polar state © is (F 2 + FT) 1 ' 2 = and that 



of the broken axisymmetry state (|8j) is (F 2 + F 2 ) 1 ! 2 = 

(l-gVSc 2 ) 172 - 

We study the stability of the polar state ([7]) using the 
Bogoliubov analysis. Substituting 
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: <fc ' r a±i,k(*), (11) 



into the GP equation (j6j), where V is the volume of the 
system, and keeping the first order terms in a±i } fc, we 
obtain 



da±i,k(t) 



dt 

where Ek 



■ cino)a±i i k(t) + cin a 



ft k 2 /(2M). The solution is given by 
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When g > g c , is real for all fc, and Eq. ([T3]) is an 
oscillating function. In this case, the polar state © is 
stable against small deviations. When q < q c , E^ is 
imaginary for < < q c — q. The modes with imaginary 
Ek exponentially grow, which make the polar state flJJ) 
dynamically unstable. 

If the initial state is prepared in the stable polar state 
with q > q c and q is decreased to g < q c , the system 
becomes dynamically unstable and the transverse mag- 
netization F± emerges 16,17 . Using Eq. ([TT]) with Eq. (fl3|h 
the correlation function of the transverse magnetization 
F± = F x ± ii 7 ^ is calculated to be 



(F + (r,i)F_(r',0) 
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where (• • • ) indicates the average with respect to dif- 
ferent initial values a±i,±fc(0). They include quantum 
and thermal fluctuations, residual atoms in the m = ±1 
states, and other experimental noises, and therefore we 
assume that a±i,±fc(0) are independent complex random 
numbers. For the dynamically unstable modes, Eq. JTSJ) 
contains the exponentially growing factor exp(2\Ek\t/h), 
which has a sharp peak at the most unstable wave num- 
ber fc mu . We can thus approximate Eq. ([T5]) as 21 
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where Ak = k — & mu , and r and £ corr are defined by 

2 -^= t -(l- 1 -e cO rAk^+0(Ak% (17) 

For q c /2 < q < q c , the most unstable wave number is 
kmu = with 



Co 
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In this case, the integral in Eq. ([T6]) for a two-dimensional 
system can be performed to yield 
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For q < <7c/2, the most unstable wave number is 
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The integral in Eq. (jl6|) for a two-dimensional system 
becomes 

{F+(r,t)F_(r',t)) oc J kJ (k\r - r '\)e-&£°" Ak * dk, 

(24) 

where Jo is the Bessel function. For k mu \r — r'\ ^> 1, this 
expression can be evaluated to be 



(F + (r,t)F_(r',t)) 

~ cos(/c mu r — 7r/4) exp 
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III. NUMERICAL RESULTS 



We restrict ourselves to two-dimensional (2D) systems 
confined in a harmonic potential Vtrap = Muo 2 (x 2 + y 2 )/2 
with uj/(2tt) = 2 Hz. When the system is tightly con- 
fined in the z direction and the thickness of the cloud 
is smaller than the spin healing length (typically a few 
micrometers), the spin dynamics are effectively 2D. We 
assume that the thickness in the z direction is ^ 1 /im, 
and use 2D interaction coefficients as c 2D = Cj/(l/im). 

We numerically solve the 2D GP equation using the 
pseudospectral method 28 . The initial state is the ground 
state of Eq. (|3]) with ip±i = 0, which is obtained by 



the imaginary time propagation method. We add small 
random noises to the initial state of ijj±i to trigger the 
magnetization. We take a sufficiently large space so that 
the boundary condition does not affect the results. 

We define the transverse and longitudinal autocorrela- 
tion functions as 



I\F + \ 2 dr ^ 
J p 2 dr ' 



J p 2 dr ' 



(26) 



We also define the transverse autocorrelation function 
along a circle with radius r as 



G T (r) 



J^\F + \ 2 (r,0)d0 
J^pHr,0)d0 ' 



(27) 



The transverse spin winding number along a circle with 
radius r is defined as 



w(r) = 



) = -/ 



d_ 
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A. Sudden quench 

We first investigate the magnetization dynamics for 
sudden quench of the quadratic Zeeman energy to q = 0. 
This corresponds to the situation in which the stable po- 
lar state is prepared at sufficiently large and the mag- 
netic field is suddenly switched off at t = 0. Figure QJa) 
shows the time evolution of the autocorrelation functions 
Gt(£) and Gl(£). The transverse magnetization starts 
to grow at t ^ 100 ms and the longitudinal magnetiza- 
tion follows. Figures [1Jb)-[3Je) show the profiles of the 
transverse magnetization. The transverse magnetization 
emerges around the center and grows outward. This is 
because the growth time in Eq. ([22]) is inversely propor- 
tional to the atomic density and the magnetization grows 
fast at which the density is large. The total density dis- 
tribution p(r) is almost unchanged during the time evo- 
lution, since Co is much larger than c\. 

Many spin vortices can be seen in Figs. ffl^-QJe) (the 
holes in the \F+\ profiles, around which argF + rotate 
by ±27r). In terms of the spin components in Eq. J5J), f3 
changes by ±27r around the vortex core, which is occupied 
by the m = component. Such a spin vortex is called a 
polar-core vortex. The spin winding number w(r) defined 
in Eq. ([28]) represents the difference between the numbers 
of polar-core vortices with opposite circulations within 
the radius r. 

We note that the spin vortices are produced by two 
distinct mechanisms in Fig. [1] with q = 0: the KZ mecha- 
nism and the spin conservation dynamic o 21 i 29 . Since the 
spin correlation function in Eq. ([25]) has a finite correla- 
tion length £com the directions of magnetization at r and 
r' are independent for \r — r'\ ^> £ C om giving rise to the 
KZ mechanism. On the other hand, when q = 0, the to- 
tal magnetization J Fdr must be conserved at zero, since 
Eq. (|3|) is invariant with respect to spin rotation in the 
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FIG. 1: (a) Time evolutions of the autocorrelation func- 
tions Gt(JL) and Gx(t) for sudden quench to q = 0. (b)-(e) 
Snapshots of the transverse magnetization \F+(r, t)\ (left pan- 
els) and its direction argF + (r,t) (right panels). The unit of 
|F+(r,t)| is 1.2 x 10 14 m" 2 . The field of view of each panel 
is 400 x 400 /xm. The number of atoms is iV = 10 7 . See the 
supplementary data file for the movie showing the dynamics. 



rotating frame ip±i — >• e^ ipt / n ip±i. In the present case, 
however, the conservation law is more strict because of 
the finite spin correlation length. Since the spin direc- 
tions at r and r' are independent for \r — r'\ ^> £ C om 
not only the total magnetization but also the local mag- 
netization J* local Fdr integrated over the size of ~ £ corr 
must be conserved in each spatial region. The magneti- 
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FIG. 2: Variance of the winding number (w 2 (r)) along the 
circumference of a circle of radius r for sudden quench, where 
the parameters are the same as those in Fig. [T] The data for 
each r is taken when Gt(v) exceeds 0.1. The average (• • •) 
is taken over 400 runs of simulations for the different initial 
states produced by random numbers. The dashed curve is the 
least square fit of Eq. (|29|) . 



zation thus occurs in such a way that the local magnetiza- 
tion is conserved at zero, i.e., spin textures are formed 30 . 
Among various spin textures, the polar-core vortices are 
favorable, since the excess energy at the defect can be 
minimized 31 . This is the second mechanism of the spin 
vortex formation in Fig. [TJ Thus, to see the effect of 
the KZ mechanism, we must take spatial region much 
larger than £ C orr- The correlation length is £ cor r — 10 
fim around the center of the trap and £ corr ~ 20 fim at 
r = 200 /im for Fig. [lj 

We consider the r-dependence of the spin winding 
number w(r). According to the KZ theory, the number of 
domains along the circle of radius r is ~ r/£ corr and hence 
w 2 (r) ~ r/£ corr . Substituting q = and q c = 2|ci|nxF(^) 
into £ corr in Eq. ([23]) , where nxF(j) oc R\ ¥ — r 2 with Rtf 



being the Thomas-Fermi radius, we obtain 



w 2 (r) oc r^. 



R 



TF 



(29) 



To compare Eq. (|29|) with the numerical simulation, we 
perform many runs of time evolution with different initial 
random noises, and take the average of w 2 (r) with respect 
to the runs, which is shown in Fig. [21 Since the time at 
which the magnetization emerges depends on r, each w(r) 
is calculated when Gt(t) in Eq. ([27]) exceeds a certain 
value (0.1 in Fig. [2]). The numerical result and Eq. ([29]) 
(circles and dashed curve in Fig. [2j respectively) are in 
good agreement, where the fitting parameter is only the 
proportionality coefficient in Eq. (|2!T 



B. Gradual quench 

We next consider the case of gradual quench of the 
magnetic field. The quadratic Zeeman energy q is linearly 
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decreased in the time scale tq as 

q(t) = qo(l-t/T Q ), 



(30) 



for < t < tq and q(t) = for t > tq. As seen in the 
previous subsection, the critical value q c (r) for magne- 
tization depends on the position r, and qo is chosen to 
be the maximum of q c (r). We define the time T(r) at 
which q reaches a local critical value as q(T(r)) = q c (r). 
The magnetization at position r is expected to emerge at 
t = T(r) + At(r) satisfying^ 



At(r) - r(r,t). 



(31) 



Using Eqs. CEJ) and §T}) with q At(r)/[q c (r)T Q ] <C 1, we 
obtain 



At(r) 



q c (r)qo 



1/3 



.1/3 



(32) 



Substituting this time into Eq. ([T9]) . we obtain the tq- 
dependence of the correlation length as 



<fcorr (f) 



Qc(r)qo 
h 2 



1/6 



1/3 



The winding number thus obeys 

2/ \ -1/3 

w (r) ex r Q . 



(33) 



(34) 



Figure [3] (a) shows (w 2 (r)} obtained by numerical sim- 
ulations of the GP equation (|6j). The variance of the 
winding number (w 2 (r)) is roughly proportional to Tq 1 ^ 3 
for 100 ms < tq < 1 s, which agrees with the above theo- 
retical argument ([34]) . For tq < 100 ms, the assumption 
of t/TQ <C 1 is violated; (w 2 (r)} approaches the values 
for sudden quench in the limit of tq — » 0. For tq > 1 

s, (w 2 (r)} significantly deviates from Tq 1 ^ 3 and steeply 
drops in Fig. [3] (a). To understand this behavior, we com- 
pare the dynamics of |F + (r)| for tq = 1 s and tq = 3 s 
shown in Figs. [3fb) and[3jc). When tq = 1 s, new spin 
vortices are produced one after another as the magneti- 
zation grows outward. When tq = 3 s, by contrast, the 
spin vortex created around the center is pushed outward, 
and no new spin vortices are created at the front of the 
magnetization growth. 

In the dynamics in Figs. [3fb) and[3^c), there are two 
characteristic velocities: the velocity v m at which the 
magnetization front spreads out and the sound velocity 
v s of the spin wave. The former is roughly obtained from 



q(t) = q c (r M (t)), 



(35) 



where tm(^) is the radius of the magnetization front. Us- 
ing the Thomas-Fermi density distribution, the right- 
hand side is q c (r M (t)) = 2|ci|n T F(rM CO) - qo[l - 



(t)/R\ ¥ \. The velocity v m is thus given by 



dr M (t) 
dt 
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(36) 
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FIG. 3: (a) Variance of the winding number (w 2 (r)) along 
the circumference of a circle of radius r = 200 fim (circles), 
r = 250 /im (squares), and r = 300 /im (triangles) for grad- 
ual quench given by Eq. (|3Q|) . The data for each r is taken 
when Gt^w) exceeds 0.1. The average (• • • ) is taken over 100 
runs of simulations for the different initial states produced by 

— 1/3 

random numbers. The dashed line is proportional to t q . 
The number of atoms is N = 10 8 . (b), (c) Snapshots of 
the transverse magnetization \F+(r,t)\ for tq = 1 s and 3 s. 
The arrows trace the vortex motion. The unit of \F+(r, t)\ is 
3.4 x 10 14 m" 2 . The field of view of each panel is 300 x 300 
/im. See the supplementary data files for the movies showing 
the dynamics. 



The transverse spin wave for the broken axisymmetry 
state ((HJ) is a phonon-like mode in the limit of k —> 0, 
whose velocity is given b y 14 i 27 



2M' 



(37) 



where q(t) in Eq. ([30]) should be used on the right-hand 
side. If v m is always faster than v 8 , the region that is go- 
ing to magnetize is causally disconnected with the mag- 
netized region, and therefore the KZ mechanism works. 
It follows from Eqs. ([36]) and ([37]) that this condition is 
satisfied for 



/2M o 

V qo 



(38) 



For the parameters in Fig. [3j the right-hand side of this 
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the potential Vtrap(^) has a minimum at r ^ 250 jam and 
the atomic density becomes maximal around this radius. 
As a result, the magnetization starts from the annulus 
around r c± 250 jam. Therefore, magnetic domains on 
this radius are always causally disconnected, and the KZ 
mechanism is expected to work even for large tq. Fig- 
ure 0] shows the results of the numerical simulations. The 
density at r ^ 250 jam is almost the same as the density 
at the same radius in the system of Fig.[3j and (w 2 (r)} is 
similar to the corresponding data in Fig. [3fa) (squares) 
for tq < 1 s. However, (w 2 (r)} in Fig. [4] obeys the KZ 

scaling Tq 1 ^ 3 even for tq > 1 s, as expected. 

IV. CONCLUSIONS 



FIG. 4: Variance of the winding number (w 2 (r)) along the 
circumference of a circle of radius r = 250 jam. (dashed circle) 
for gradual quench given by Eq. (|30|) . The potential has a 
form of Eq. (|39|) . The data for each r is taken when Gt(v w ) 
exceeds 0.1. The average (• • • ) is taken over 100 runs of sim- 
ulations for the different initial states produced by random 
numbers. The dashed line is least square fit of the plots by 

2 /3 

using a function proportional to Tq 1 . The number of atoms 
is N = 10 8 . The inset shows the initial density profile with the 
unit of density 2.4 x 10 14 m~ 2 and the field of view 900 x 900 
/am. 

inequality is ~ 1.7 s, which agrees well with the time at 
which the plots in Fig. [3fa) deviates from t q 1 . 

C. Gradual quench with a plug potential 



We have investigated the spin vortex formation due 
to the KZ mechanism in a quenched ferromagnetic BEC 
confined in a trapping potential. Since the atomic density 
is inhomogeneous in a harmonic trap, the spin correlation 
length depends on the radius r. In fact, the numerical 
simulations showed that the spin winding number de- 
pends on r, which was in good agreement with the theo- 
retical prediction (Fig. [2]). When the quadratic Zeeman 
energy q is gradually quenched with the time scale tq , the 
magnetized region gradually expands from the center to 
the periphery of the atomic cloud. If the expansion veloc- 
ity is much faster than the spin wave velocity, the system 
exhibits the KZ scaling law, and if the former is slower 
than the latter, the KZ scenario breaks down (Fig. [3j). 
When a plug potential is added to the harmonic trap, 
the geometry of the system is changed and the KZ power 
law can be observed over a wide range of tq (Fig. [4]). 



In the previous subsection, we showed that the KZ sce- 
nario breaks down when the magnetized region expands 
slowly. This is because the region that is going to mag- 
netize is causally connected to the initially magnetized 
region at the trap center. To eliminate this effect, we 
remove atoms around the trap center by adding a plug 
potential as 

V trap (r) = IjlfwV + Ae- r *' d \ (39) 

where the values of the parameters are chosen to be 
A = 1500fkj and d = 222 jam. For these parameters, 
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